home *** CD-ROM | disk | FTP | other *** search
/ Chip 2007 January, February, March & April / Chip-Cover-CD-2007-02.iso / Pakiet bezpieczenstwa / mini Pentoo LiveCD 2006.1 / mpentoo-2006.1.iso / livecd.squashfs / usr / lib / python2.4 / site-packages / Numeric / LinearAlgebra.pyc (.txt) < prev    next >
Encoding:
Python Compiled Bytecode  |  2006-03-29  |  15.1 KB  |  513 lines
  1. # Source Generated with Decompyle++
  2. # File: in.pyc (Python 2.4)
  3.  
  4. import Numeric
  5. import copy
  6. import lapack_lite
  7. import math
  8. import MLab
  9. import multiarray
  10.  
  11. class LinAlgError(Exception):
  12.     pass
  13.  
  14. _lapack_type = {
  15.     'f': 0,
  16.     'd': 1,
  17.     'F': 2,
  18.     'D': 3 }
  19. _lapack_letter = [
  20.     's',
  21.     'd',
  22.     'c',
  23.     'z']
  24. _array_kind = {
  25.     'i': 0,
  26.     'l': 0,
  27.     'f': 0,
  28.     'd': 0,
  29.     'F': 1,
  30.     'D': 1 }
  31. _array_precision = {
  32.     'i': 1,
  33.     'l': 1,
  34.     'f': 0,
  35.     'd': 1,
  36.     'F': 0,
  37.     'D': 1 }
  38. _array_type = [
  39.     [
  40.         'f',
  41.         'd'],
  42.     [
  43.         'F',
  44.         'D']]
  45.  
  46. def _commonType(*arrays):
  47.     kind = 0
  48.     precision = 1
  49.     for a in arrays:
  50.         t = a.typecode()
  51.         kind = max(kind, _array_kind[t])
  52.         precision = max(precision, _array_precision[t])
  53.     
  54.     return _array_type[kind][precision]
  55.  
  56.  
  57. def _castCopyAndTranspose(type, *arrays):
  58.     cast_arrays = ()
  59.     for a in arrays:
  60.         if a.typecode() == type:
  61.             cast_arrays = cast_arrays + (copy.copy(Numeric.transpose(a)),)
  62.             continue
  63.         cast_arrays = cast_arrays + (copy.copy(Numeric.transpose(a).astype(type)),)
  64.     
  65.     if len(cast_arrays) == 1:
  66.         return cast_arrays[0]
  67.     else:
  68.         return cast_arrays
  69.  
  70. _fastCT = multiarray._fastCopyAndTranspose
  71.  
  72. def _fastCopyAndTranspose(type, *arrays):
  73.     cast_arrays = ()
  74.     for a in arrays:
  75.         if a.typecode() == type:
  76.             cast_arrays = cast_arrays + (_fastCT(a),)
  77.             continue
  78.         cast_arrays = cast_arrays + (_fastCT(a.astype(type)),)
  79.     
  80.     if len(cast_arrays) == 1:
  81.         return cast_arrays[0]
  82.     else:
  83.         return cast_arrays
  84.  
  85.  
  86. def _assertRank2(*arrays):
  87.     for a in arrays:
  88.         if len(a.shape) != 2:
  89.             raise LinAlgError, 'Array must be two-dimensional'
  90.             continue
  91.     
  92.  
  93.  
  94. def _assertSquareness(*arrays):
  95.     for a in arrays:
  96.         if max(a.shape) != min(a.shape):
  97.             raise LinAlgError, 'Array must be square'
  98.             continue
  99.     
  100.  
  101.  
  102. def solve_linear_equations(a, b):
  103.     one_eq = len(b.shape) == 1
  104.     if one_eq:
  105.         b = b[(:, Numeric.NewAxis)]
  106.     
  107.     _assertRank2(a, b)
  108.     _assertSquareness(a)
  109.     n_eq = a.shape[0]
  110.     n_rhs = b.shape[1]
  111.     if n_eq != b.shape[0]:
  112.         raise LinAlgError, 'Incompatible dimensions'
  113.     
  114.     t = _commonType(a, b)
  115.     if _array_kind[t] == 1:
  116.         lapack_routine = lapack_lite.zgesv
  117.     else:
  118.         lapack_routine = lapack_lite.dgesv
  119.     (a, b) = _fastCopyAndTranspose(t, a, b)
  120.     pivots = Numeric.zeros(n_eq, 'i')
  121.     results = lapack_routine(n_eq, n_rhs, a, n_eq, pivots, b, n_eq, 0)
  122.     if results['info'] > 0:
  123.         raise LinAlgError, 'Singular matrix'
  124.     
  125.     if one_eq:
  126.         return Numeric.ravel(b)
  127.     else:
  128.         return multiarray.transpose(b)
  129.  
  130.  
  131. def inverse(a):
  132.     return solve_linear_equations(a, Numeric.identity(a.shape[0]))
  133.  
  134.  
  135. def cholesky_decomposition(a):
  136.     _assertRank2(a)
  137.     _assertSquareness(a)
  138.     t = _commonType(a)
  139.     a = _castCopyAndTranspose(t, a)
  140.     m = a.shape[0]
  141.     n = a.shape[1]
  142.     if _array_kind[t] == 1:
  143.         lapack_routine = lapack_lite.zpotrf
  144.     else:
  145.         lapack_routine = lapack_lite.dpotrf
  146.     results = lapack_routine('L', n, a, m, 0)
  147.     if results['info'] > 0:
  148.         raise LinAlgError, 'Matrix is not positive definite - Cholesky decomposition cannot be computed'
  149.     
  150.     return copy.copy(Numeric.transpose(MLab.triu(a, k = 0)))
  151.  
  152.  
  153. def eigenvalues(a):
  154.     _assertRank2(a)
  155.     _assertSquareness(a)
  156.     t = _commonType(a)
  157.     real_t = _array_type[0][_array_precision[t]]
  158.     a = _fastCopyAndTranspose(t, a)
  159.     n = a.shape[0]
  160.     dummy = Numeric.zeros((1,), t)
  161.     if _array_kind[t] == 1:
  162.         lapack_routine = lapack_lite.zgeev
  163.         w = Numeric.zeros((n,), t)
  164.         rwork = Numeric.zeros((n,), real_t)
  165.         lwork = 1
  166.         work = Numeric.zeros((lwork,), t)
  167.         results = lapack_routine('N', 'N', n, a, n, w, dummy, 1, dummy, 1, work, -1, rwork, 0)
  168.         lwork = int(abs(work[0]))
  169.         work = Numeric.zeros((lwork,), t)
  170.         results = lapack_routine('N', 'N', n, a, n, w, dummy, 1, dummy, 1, work, lwork, rwork, 0)
  171.     else:
  172.         lapack_routine = lapack_lite.dgeev
  173.         wr = Numeric.zeros((n,), t)
  174.         wi = Numeric.zeros((n,), t)
  175.         lwork = 1
  176.         work = Numeric.zeros((lwork,), t)
  177.         results = lapack_routine('N', 'N', n, a, n, wr, wi, dummy, 1, dummy, 1, work, -1, 0)
  178.         lwork = int(work[0])
  179.         work = Numeric.zeros((lwork,), t)
  180.         results = lapack_routine('N', 'N', n, a, n, wr, wi, dummy, 1, dummy, 1, work, lwork, 0)
  181.         if Numeric.logical_and.reduce(Numeric.equal(wi, 0.0)):
  182.             w = wr
  183.         else:
  184.             w = wr + (0.0+1.0j) * wi
  185.     if results['info'] > 0:
  186.         raise LinAlgError, 'Eigenvalues did not converge'
  187.     
  188.     return w
  189.  
  190.  
  191. def Heigenvalues(a, UPLO = 'L'):
  192.     _assertRank2(a)
  193.     _assertSquareness(a)
  194.     t = _commonType(a)
  195.     real_t = _array_type[0][_array_precision[t]]
  196.     a = _castCopyAndTranspose(t, a)
  197.     n = a.shape[0]
  198.     liwork = 5 * n + 3
  199.     iwork = Numeric.zeros((liwork,), 'i')
  200.     if _array_kind[t] == 1:
  201.         lapack_routine = lapack_lite.zheevd
  202.         w = Numeric.zeros((n,), real_t)
  203.         lwork = 1
  204.         work = Numeric.zeros((lwork,), t)
  205.         lrwork = 1
  206.         rwork = Numeric.zeros((lrwork,), real_t)
  207.         results = lapack_routine('N', UPLO, n, a, n, w, work, -1, rwork, -1, iwork, liwork, 0)
  208.         lwork = int(abs(work[0]))
  209.         work = Numeric.zeros((lwork,), t)
  210.         lrwork = int(rwork[0])
  211.         rwork = Numeric.zeros((lrwork,), real_t)
  212.         results = lapack_routine('N', UPLO, n, a, n, w, work, lwork, rwork, lrwork, iwork, liwork, 0)
  213.     else:
  214.         lapack_routine = lapack_lite.dsyevd
  215.         w = Numeric.zeros((n,), t)
  216.         lwork = 1
  217.         work = Numeric.zeros((lwork,), t)
  218.         results = lapack_routine('N', UPLO, n, a, n, w, work, -1, iwork, liwork, 0)
  219.         lwork = int(work[0])
  220.         work = Numeric.zeros((lwork,), t)
  221.         results = lapack_routine('N', UPLO, n, a, n, w, work, lwork, iwork, liwork, 0)
  222.     if results['info'] > 0:
  223.         raise LinAlgError, 'Eigenvalues did not converge'
  224.     
  225.     return w
  226.  
  227.  
  228. def eigenvectors(a):
  229.     '''eigenvectors(a) returns u,v  where u is the eigenvalues and
  230. v is a matrix of eigenvectors with vector v[i] corresponds to 
  231. eigenvalue u[i].  Satisfies the equation dot(a, v[i]) = u[i]*v[i]
  232. '''
  233.     _assertRank2(a)
  234.     _assertSquareness(a)
  235.     t = _commonType(a)
  236.     real_t = _array_type[0][_array_precision[t]]
  237.     a = _fastCopyAndTranspose(t, a)
  238.     n = a.shape[0]
  239.     dummy = Numeric.zeros((1,), t)
  240.     if _array_kind[t] == 1:
  241.         lapack_routine = lapack_lite.zgeev
  242.         w = Numeric.zeros((n,), t)
  243.         v = Numeric.zeros((n, n), t)
  244.         lwork = 1
  245.         work = Numeric.zeros((lwork,), t)
  246.         rwork = Numeric.zeros((2 * n,), real_t)
  247.         results = lapack_routine('N', 'V', n, a, n, w, dummy, 1, v, n, work, -1, rwork, 0)
  248.         lwork = int(abs(work[0]))
  249.         work = Numeric.zeros((lwork,), t)
  250.         results = lapack_routine('N', 'V', n, a, n, w, dummy, 1, v, n, work, lwork, rwork, 0)
  251.     else:
  252.         lapack_routine = lapack_lite.dgeev
  253.         wr = Numeric.zeros((n,), t)
  254.         wi = Numeric.zeros((n,), t)
  255.         vr = Numeric.zeros((n, n), t)
  256.         lwork = 1
  257.         work = Numeric.zeros((lwork,), t)
  258.         results = lapack_routine('N', 'V', n, a, n, wr, wi, dummy, 1, vr, n, work, -1, 0)
  259.         lwork = int(work[0])
  260.         work = Numeric.zeros((lwork,), t)
  261.         results = lapack_routine('N', 'V', n, a, n, wr, wi, dummy, 1, vr, n, work, lwork, 0)
  262.         if Numeric.logical_and.reduce(Numeric.equal(wi, 0.0)):
  263.             w = wr
  264.             v = vr
  265.         else:
  266.             w = wr + (0.0+1.0j) * wi
  267.             v = Numeric.array(vr, Numeric.Complex)
  268.             ind = Numeric.nonzero(Numeric.equal(Numeric.equal(wi, 0.0), 0))
  269.             for i in range(len(ind) / 2):
  270.                 v[ind[2 * i]] = vr[ind[2 * i]] + (0.0+1.0j) * vr[ind[2 * i + 1]]
  271.                 v[ind[2 * i + 1]] = vr[ind[2 * i]] - (0.0+1.0j) * vr[ind[2 * i + 1]]
  272.             
  273.     if results['info'] > 0:
  274.         raise LinAlgError, 'Eigenvalues did not converge'
  275.     
  276.     return (w, v)
  277.  
  278.  
  279. def Heigenvectors(a, UPLO = 'L'):
  280.     _assertRank2(a)
  281.     _assertSquareness(a)
  282.     t = _commonType(a)
  283.     real_t = _array_type[0][_array_precision[t]]
  284.     a = _castCopyAndTranspose(t, a)
  285.     n = a.shape[0]
  286.     liwork = 5 * n + 3
  287.     iwork = Numeric.zeros((liwork,), 'i')
  288.     if _array_kind[t] == 1:
  289.         lapack_routine = lapack_lite.zheevd
  290.         w = Numeric.zeros((n,), real_t)
  291.         lwork = 1
  292.         work = Numeric.zeros((lwork,), t)
  293.         lrwork = 1
  294.         rwork = Numeric.zeros((lrwork,), real_t)
  295.         results = lapack_routine('V', UPLO, n, a, n, w, work, -1, rwork, -1, iwork, liwork, 0)
  296.         lwork = int(abs(work[0]))
  297.         work = Numeric.zeros((lwork,), t)
  298.         lrwork = int(rwork[0])
  299.         rwork = Numeric.zeros((lrwork,), real_t)
  300.         results = lapack_routine('V', UPLO, n, a, n, w, work, lwork, rwork, lrwork, iwork, liwork, 0)
  301.     else:
  302.         lapack_routine = lapack_lite.dsyevd
  303.         w = Numeric.zeros((n,), t)
  304.         lwork = 1
  305.         work = Numeric.zeros((lwork,), t)
  306.         results = lapack_routine('V', UPLO, n, a, n, w, work, -1, iwork, liwork, 0)
  307.         lwork = int(work[0])
  308.         work = Numeric.zeros((lwork,), t)
  309.         results = lapack_routine('V', UPLO, n, a, n, w, work, lwork, iwork, liwork, 0)
  310.     if results['info'] > 0:
  311.         raise LinAlgError, 'Eigenvalues did not converge'
  312.     
  313.     return (w, a)
  314.  
  315.  
  316. def singular_value_decomposition(a, full_matrices = 0):
  317.     _assertRank2(a)
  318.     n = a.shape[1]
  319.     m = a.shape[0]
  320.     t = _commonType(a)
  321.     real_t = _array_type[0][_array_precision[t]]
  322.     a = _fastCopyAndTranspose(t, a)
  323.     if full_matrices:
  324.         nu = m
  325.         nvt = n
  326.         option = 'A'
  327.     else:
  328.         nu = min(n, m)
  329.         nvt = min(n, m)
  330.         option = 'S'
  331.     s = Numeric.zeros((min(n, m),), real_t)
  332.     u = Numeric.zeros((nu, m), t)
  333.     vt = Numeric.zeros((n, nvt), t)
  334.     iwork = Numeric.zeros((8 * min(m, n),), 'i')
  335.     if _array_kind[t] == 1:
  336.         lapack_routine = lapack_lite.zgesdd
  337.         rwork = Numeric.zeros((5 * min(m, n) * min(m, n) + 5 * min(m, n),), real_t)
  338.         lwork = 1
  339.         work = Numeric.zeros((lwork,), t)
  340.         results = lapack_routine(option, m, n, a, m, s, u, m, vt, nvt, work, -1, rwork, iwork, 0)
  341.         lwork = int(abs(work[0]))
  342.         work = Numeric.zeros((lwork,), t)
  343.         results = lapack_routine(option, m, n, a, m, s, u, m, vt, nvt, work, lwork, rwork, iwork, 0)
  344.     else:
  345.         lapack_routine = lapack_lite.dgesdd
  346.         lwork = 1
  347.         work = Numeric.zeros((lwork,), t)
  348.         results = lapack_routine(option, m, n, a, m, s, u, m, vt, nvt, work, -1, iwork, 0)
  349.         lwork = int(work[0])
  350.         work = Numeric.zeros((lwork,), t)
  351.         results = lapack_routine(option, m, n, a, m, s, u, m, vt, nvt, work, lwork, iwork, 0)
  352.     if results['info'] > 0:
  353.         raise LinAlgError, 'SVD did not converge'
  354.     
  355.     return (multiarray.transpose(u), s, multiarray.transpose(vt))
  356.  
  357.  
  358. def generalized_inverse(a, rcond = 1e-10):
  359.     a = Numeric.array(a, copy = 0)
  360.     if a.typecode() in Numeric.typecodes['Complex']:
  361.         a = Numeric.conjugate(a)
  362.     
  363.     (u, s, vt) = singular_value_decomposition(a, 0)
  364.     m = u.shape[0]
  365.     n = vt.shape[1]
  366.     cutoff = rcond * Numeric.maximum.reduce(s)
  367.     for i in range(min(n, m)):
  368.         if s[i] > cutoff:
  369.             s[i] = 1.0 / s[i]
  370.             continue
  371.         s[i] = 0.0
  372.     
  373.     return Numeric.dot(Numeric.transpose(vt), s[(:, Numeric.NewAxis)] * Numeric.transpose(u))
  374.  
  375.  
  376. def determinant(a):
  377.     _assertRank2(a)
  378.     _assertSquareness(a)
  379.     t = _commonType(a)
  380.     a = _fastCopyAndTranspose(t, a)
  381.     n = a.shape[0]
  382.     if _array_kind[t] == 1:
  383.         lapack_routine = lapack_lite.zgetrf
  384.     else:
  385.         lapack_routine = lapack_lite.dgetrf
  386.     pivots = Numeric.zeros((n,), 'i')
  387.     results = lapack_routine(n, n, a, n, pivots, 0)
  388.     sign = Numeric.add.reduce(Numeric.not_equal(pivots, Numeric.arrayrange(1, n + 1))) % 2
  389.     return (1.0 - 2.0 * sign) * Numeric.multiply.reduce(Numeric.diagonal(a))
  390.  
  391.  
  392. def linear_least_squares(a, b, rcond = 1e-10):
  393.     '''solveLinearLeastSquares(a,b) returns x,resids,rank,s 
  394. where x minimizes 2-norm(|b - Ax|) 
  395.       resids is the sum square residuals
  396.       rank is the rank of A
  397.       s is an rank of the singual values of A in desending order
  398.  
  399. If b is a matrix then x is also a matrix with corresponding columns.
  400. If the rank of A is less than the number of columns of A or greater than
  401. the numer of rows, then residuals will be returned as an empty array
  402. otherwise resids = sum((b-dot(A,x)**2).
  403. Singular values less than s[0]*rcond are treated as zero.
  404. '''
  405.     one_eq = len(b.shape) == 1
  406.     if one_eq:
  407.         b = b[(:, Numeric.NewAxis)]
  408.     
  409.     _assertRank2(a, b)
  410.     m = a.shape[0]
  411.     n = a.shape[1]
  412.     n_rhs = b.shape[1]
  413.     ldb = max(n, m)
  414.     if m != b.shape[0]:
  415.         raise LinAlgError, 'Incompatible dimensions'
  416.     
  417.     t = _commonType(a, b)
  418.     real_t = _array_type[0][_array_precision[t]]
  419.     bstar = Numeric.zeros((ldb, n_rhs), t)
  420.     bstar[(:b.shape[0], :n_rhs)] = copy.copy(b)
  421.     (a, bstar) = _castCopyAndTranspose(t, a, bstar)
  422.     s = Numeric.zeros((min(m, n),), real_t)
  423.     nlvl = max(0, int(math.log(float(min(m, n)) / 2.0)) + 1)
  424.     iwork = Numeric.zeros((3 * min(m, n) * nlvl + 11 * min(m, n),), 'i')
  425.     if _array_kind[t] == 1:
  426.         lapack_routine = lapack_lite.zgelsd
  427.         lwork = 1
  428.         rwork = Numeric.zeros((lwork,), real_t)
  429.         work = Numeric.zeros((lwork,), t)
  430.         results = lapack_routine(m, n, n_rhs, a, m, bstar, ldb, s, rcond, 0, work, -1, rwork, iwork, 0)
  431.         lwork = int(abs(work[0]))
  432.         rwork = Numeric.zeros((lwork,), real_t)
  433.         a_real = Numeric.zeros((m, n), real_t)
  434.         bstar_real = Numeric.zeros((ldb, n_rhs), real_t)
  435.         results = lapack_lite.dgelsd(m, n, n_rhs, a_real, m, bstar_real, ldb, s, rcond, 0, rwork, -1, iwork, 0)
  436.         lrwork = int(rwork[0])
  437.         work = Numeric.zeros((lwork,), t)
  438.         rwork = Numeric.zeros((lrwork,), real_t)
  439.         results = lapack_routine(m, n, n_rhs, a, m, bstar, ldb, s, rcond, 0, work, lwork, rwork, iwork, 0)
  440.     else:
  441.         lapack_routine = lapack_lite.dgelsd
  442.         lwork = 1
  443.         work = Numeric.zeros((lwork,), t)
  444.         results = lapack_routine(m, n, n_rhs, a, m, bstar, ldb, s, rcond, 0, work, -1, iwork, 0)
  445.         lwork = int(work[0])
  446.         work = Numeric.zeros((lwork,), t)
  447.         results = lapack_routine(m, n, n_rhs, a, m, bstar, ldb, s, rcond, 0, work, lwork, iwork, 0)
  448.     if results['info'] > 0:
  449.         raise LinAlgError, 'SVD did not converge in Linear Least Squares'
  450.     
  451.     resids = Numeric.array([], t)
  452.     if one_eq:
  453.         x = copy.copy(Numeric.ravel(bstar)[:n])
  454.         if results['rank'] == n and m > n:
  455.             resids = Numeric.array([
  456.                 Numeric.sum(Numeric.ravel(bstar)[n:] ** 2)])
  457.         
  458.     else:
  459.         x = copy.copy(Numeric.transpose(bstar)[(:n, :)])
  460.         if results['rank'] == n and m > n:
  461.             resids = copy.copy(Numeric.sum(Numeric.transpose(bstar)[(n:, :)] ** 2))
  462.         
  463.     return (x, resids, results['rank'], copy.copy(s[:min(n, m)]))
  464.  
  465. if __name__ == '__main__':
  466.     from Numeric import *
  467.     
  468.     def test(a, b):
  469.         print 'All numbers printed should be (almost) zero:'
  470.         x = solve_linear_equations(a, b)
  471.         check = b - matrixmultiply(a, x)
  472.         print check
  473.         a_inv = inverse(a)
  474.         check = matrixmultiply(a, a_inv) - identity(a.shape[0])
  475.         print check
  476.         ev = eigenvalues(a)
  477.         (evalues, evectors) = eigenvectors(a)
  478.         check = ev - evalues
  479.         print check
  480.         evectors = transpose(evectors)
  481.         check = matrixmultiply(a, evectors) - evectors * evalues
  482.         print check
  483.         (u, s, vt) = singular_value_decomposition(a)
  484.         check = a - Numeric.matrixmultiply(u * s, vt)
  485.         print check
  486.         a_ginv = generalized_inverse(a)
  487.         check = matrixmultiply(a, a_ginv) - identity(a.shape[0])
  488.         print check
  489.         det = determinant(a)
  490.         check = det - multiply.reduce(evalues)
  491.         print check
  492.         (x, residuals, rank, sv) = linear_least_squares(a, b)
  493.         check = b - matrixmultiply(a, x)
  494.         print check
  495.         print rank - a.shape[0]
  496.         print sv - s
  497.  
  498.     a = array([
  499.         [
  500.             1.0,
  501.             2.0],
  502.         [
  503.             3.0,
  504.             4.0]])
  505.     b = array([
  506.         2.0,
  507.         1.0])
  508.     test(a, b)
  509.     a = a + (0.0+0.0j)
  510.     b = b + (0.0+0.0j)
  511.     test(a, b)
  512.  
  513.